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Preface 



The simulation of the dynamics of classical field theories is currently gaining some attention 
from the high-energy community, mainly in the context of statistical field theory. Recent papers 
show that, in some particular but important conditions, classical field theories are a very good 
approximation to the quantum evolution of fields at finite temperature (see, for instance, |j], 
H]). Also, in the context of condensed-matter systems, the dynamics of effective classical fields 
has proven to be a very efficient tool in describing both the equilibrium and non-equilibrium 
properties of systems such as ferromagnets ||. In order to simulate these theories, they must 
be first discretized and cast on a lattice, a job far from trivial to do in a consistent manner 
||. After the discretization, the differential equations of motion transform themselves into finite 
difference equations. Before doing any sort of useful calculation in the physical context above, one 
is supposed to understand the basic foundations of the numerical method, the study of which is 
the objective of this monograph. 

Following this motivation, I will present, in an introductory way, the Finite Difference method 
for hyperbolic equations, focusing on a method which has second order precision both in time and 
space (the so-called staggered leapfrog method) and applying it to the case of the Id and 2d wave 
equation. A brief derivation of the energy and equation of motion of a wave is done before the 
numerical part in order to make the transition from the continuum to the lattice clearer. 

To illustrate the extension of the method to more complex equations, I also add dissipative 
terms of the kind — r/u into the equations. I also briefly discuss the von Neumann numerical 
stability analysis and the Courant criterion, two of the most popular in the literature. In the end 
I present some numerical results obtained with the leapfrog algorithm, illustrating the importance 
of the lattice resolution through energy plots. 

I have tried to collect, in a concise way, the main steps necessary to have a stable algorithm to 
solve wave- like equations. More sophisticated versions of these equations should be handled with 
care, and accompanied of a rigorous study of convergence and stability which could be found in 
the references cited in the end of this work. 
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Chapter 1 

The Wave Equation 



1.1 Introduction 

Partial Differential Equations (from now on simply PDEs) are divided in the literature basically 
in three kinds: parabolic, elliptic and hyperbolic (the criterion of classification of these equations 
can be found in |g, Chap. 8]). In this work we will be interested mainly on hyperbolic equations, 
of which the wave equation is the paradigm: 



1 d 2 u 

lA ii — - 

„2 Q t 2 



v 2 « = 4^ (i-i) 



where v — ? is the square of the wave velocity in the medium, which in the case of a free string 
could be determined by the tension r and the mass density per length unit A. 

From the strict numerical point of view, the distinction between these classes of PDEs isn't of 
much importance ||. There is, however, another sort of classification of PDEs which is relevant for 
numerical purposes: the initial value problems (which include the case of the hyperbolic equations) 
and the boundary condition problems (which include, for instance, parabolic equations). In this 
work we will restrict ourselves to initial value problems. See reference |(j|] for a good introduction 
to boundary conditi on p roblems. 



In the equation ( |l.l| ) we could still add a dissipative term proportional to the first power of 



the time derivative of u, i.e., 



,-j , d 2 u du 



where n is the viscosity coefficient. 

Our first step will be to derive the wave equation from a simple mechanical analysis of the free 
rope. Being this a well known problem in classical mechanics, we will go through only the main 
steps of it (for a more complete treatment of the problem of the free string, see, for example, fiL 
Chaps. 8 and 9]). Once we are done with the 1-d wave equation, we will proceed further to the 
2-d case, which isn't as abundant in the literature as the 1-d case. 

1.2 Waves in 1-dimension (the free string) 

1.2.1 Equation of Motion 

Figure [l]l] gives us an idea of a mass element dm with linear dimension dx subject to tension 
forces. We are interested on the vertical displacement of this mass element, so, for this direction, 
we could write the resulting force: 

dF u = f ■u\ x+dx - t-u\ x , (1.3) 
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x x+dx 

Figure 1.1: Representing the tension forces acting on an infinitesimal element of the rope. 



where f is the tension and the unit vector u refers to the vertical direction. Within the domain of 
smooth deformations of the string (i.e., small (3), we could write: 

t)ii 



f ■ u — t sin (3 m t tan (3 = r 



dx 



(1.4) 



We notice now that ( [T^ ) could be written as: 

T iAx+dx 



, lx+ax -r u \ _ dr^ _ ^ f du 

ax ox ox \ ox . 



We will now restrict ourselves to the case of constant tensions along the rope, so that: 

jv Q2u j 

dt u = T—^dx 

ox 1 
Equating this with Newton's second law 

dF u = dm— = X-^jdx, 
where A is the linear mass density, we obtain then the wave equation for a free string: 



(1.5) 



dx 2 dt 2 



(1.6) 



1.2.2 Energy 



The kinetic energy could be evaluated in a straightforward manner, integrating the kinetic energy 
term for a representative mass element: 

1 



dT = -dm ■ v,,, 
2 u 



with dm = X ■ dx and v u = S, in other words 



da 



dt ' 



T = I dT = ± f (^\ ,1, 



(1.7) 



The potential energy can be obtained by calculating the work necessary to bring the string 
from a "trivial" configuration u(x, 0) = to the configuration at which we want to evaluate the 
potential energy u(x,t). We will fix the boundary conditions u(0,t) = u(l,t) = (string tied at 
the ends) and as a consequence of this, d t u(0,t) — d t u(l,t) = 0. The potential energy relative to 
the work necessary to change of 8u the configuration of an element of the string in an interval of 
time dt is: 

SV = -dF u ■ Su = -dF u ■(foUt 
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Therefore, the potential energy of the whole string in this same interval of time is: 

•IV = -'II- I dF u (^ 



Substituting (1.5) in the latter we get: 



„,__«. / , T (gV£ , i* 



Of 



We are however interested on V [u(x, t)], so integrating in time and using the boundary condi- 
tions above, we have: 



V = 



dV = -t dt dx 
Jo Jo 

rt , f du du ' 
t I dtl — — 



dt ) \ dx 2 



dt 



dt dx 
du d 2 u 



u 



du d' 
(l jo dx dxdt 

dx = 



dx 



Jo dx dxdt 
r * 1 d f l fdu\ 2 , 
' d Wo fej dX 
1 'du^ 2 ' 



(I 



II 



— I dx 

dx 



v -UM}- 



With ( pT^ ) and (L8) we have finally the total energy of the rope: 



A f l fdu^ 2 



E =2jA^) dx+ U!(M) 2dx 



(U 



(1.9) 



In our applications we will take X — t such that v — 1 and (|l.9|) assumes the simple form: 



E 



2 7o 






du 

dx 



dx 



(1.10) 



This energy equation will be very useful to test our algorithms through an analysis of conser- 
vation (or dissipation) during the dynamical evolution of the system. 

1.3 Waves in 2-dimensions (the free membrane) 

1.3.1 Equation of Motion 

In its two dimensional version, the wave equation could be describing a membrane, a liquid surface, 
or some "coarse-grained" field in the surface physics, to cite a few. In the case of the membrane 
or other elastic surface, the oscillations are also constrained to be small (analogously to the 1-d 
string). 

We notice now that the additional dimension forces us to define the tension "per unit length": 



f 



(1.11) 
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dm 





dy 



/ dx / 
x x+dx 



Figure 1.2: Representing a mass element dm of dimensions dxdy. The element is subject to tension 
forces on each side (analogous to the borders of the mass element in the 1-d case), being these 
forces orthogonal to the axes of the sides. Only the variation with respect to x of u is drawn (i.e., 
u and u + du in the figure are displacements of u(x, y) keeping y constant and varying x). 



This force per unit length could be understood with a simple example: stretch a tape of width 
I from its extremities with force r. We can't ask the force on a point of the tape, but only on some 
element of some definite length and width (of course, this element could be differential, playing 
the same role of a linear differential element in the case of the string) . Formula ( |l.ll| ) , times the 
length of the element, then gives you the resulting force (tension) on the element. In this way, we 
extend the equation (|1.3|) to two dimensions: 



dF u = 



JX 



n dy 

x-\-dx,y 



Jx 



dy 

x,y 



!, 



ii dx 

x,y+dy 



fy 



a dx 
x,y 



(1.12) 



where f x ■ u dy = f x . u \ dy is the x component of the tension in the direction u acting on 

x,y 

the side defined by the points (x,y) e (x,y + dy), and so on. We are going to suppose that the 
forces on the sides of the elements are orthogonal to their axes, which is the same as decomposing 
the tension force on dm into four components, one for each side (notice, however, that we have 
effectively only two resulting components, to wit x and y). Doing this we won't need to emphasize 

becoming simply / u | and so on. Nevertheless 



%,y 



the tension components along x or y, f x 

it is still important, for what we said above, to know what side we are talking about. So, in the 
regime of small vibrations (i.e., small angles of deformation), we could find f u analogously to the 
string case 



fudy = f-g^dy 

f u dx = f-^-dx, 
dy 



(1.13) 
(1.14) 



where f u dy and f u dx are the tensions in the direction u on a side of length dy and dx along y 
and x, respectively. We emphasize that, with this notation plus the knowledge of the point where 
we are going to evaluate the derivatives, we have a complete specification of the side on which the 
tension acts 1 . With all this in hands, Eq. ( 1.12J ) becomes: 



dF u = 



J u \x, 



y+dy 



dx ~ U\ x , y dx 



Ju\ x +dx,y dy Ju\ Xt ydy 

With calculations analogous to those of the previous section, we have 

J u \x-\-dx.y J u \x^y j , J u \x,y-\-dy J u \x,y 



dF„ 



dx 



dy 



dydx 



(1.15) 



1 Indeed, once specified the beginning of the side with the pair (x, y), specifying the length with dx or dy furnishes 
us with the direction of the side in question. This is sufficient to localize it, since the z coordinate is unambiguously 
determined via z = u(x,y). 



CHAPTER 1. THE WAVE EQUATION 



— — dxdy H — -—dxdy 
ox ay 



f d^ dxdy + f d^ dxdy 



With Newton's second law we obtain: 



(1.16) 
(1.17) 



f-d^ dxd y + f-Q-2 dxd y 

t^ U A A j. ^A A 

J ^15 dxdy + J -q^ dxdy 



= dm 



dt 2 



dx 2 



f d 2 u 
(dx 2 



a u 

dy 2 



A A ^ U 

adxdy-^j 



dt 2 



(1.18) 



d 2 u d 2 u 1 d 2 u 
dx 2 dy 2 v 2 dt 2 ' 

which is the desired wave equation for two dimensions, with v 2 = £ and a the surface mass density. 

1.3.2 Energy 

The derivation of the total energy is done in the same manner as the Id case. We will consider 
a surface z = u(x,y,i) with support of dimension I x I, subject to the boundary conditions 

U \boundary = u (^,y,t) = u(l,y,t) = U (x,0,t) = u(x,l,t) = and u\boundary = u(°> Vi *) = 

u(l, y, t) — ii(x, 0, t) — u(x, I, t) = 0, where u = du/dt. Let us begin with the kinetic term: 



dT 

T 
where a is the surface mass density. 



— dm ■ u = —u dxdy 
2 2 y 



l A 

o Jo 



■jjj ) ' ,X( 'H- 



(1.19) 

(1.20) 



'"■<t> = 



The potential energy is obtained in an analogous way to the Section (1.2) 
dV = -dt ■ 
= -dt- 
V = -f 



l A 



( d 2 u d 2 u\ ( ' du\ 

f , f /"' , f l , (d 2 udu\ r l , r l , (d 2 udu 

o dt \jo dy Jo dX {^) + Jo dX Jo dy {-o7m 



- nAfAUm 2 - 



/ 


I'A 


[Jo Jo 


[( 


du\ 2 

dx) 


( dv 


r 


da 


/ 


ff 

Jo Jo 


( du\ ( du\ 
\dx) \dy) 


dxdy 


t 






V 


Aj 


f l 


~(du s 
\dx, 


>+( 


dy, 


)] 



dx- — 
o 2dt J0 



dxdy > = 



%)■"» = 



tt- + -7T~ dxdy 



(1.21) 
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Then the total energy for our usual condition / = a 

r l rl 



E 



1 

2 Jo Jo 



v 
dy, 



1 is: 



i~\ \ 2, / r\ \ 2 / r\ \ 2 

ou\ I ou\ / ou 



dxdy 



v dx J \dy J \ dt t 
or, in a more general way (we are going to take as granted the result for more than two dimensions), 

~1,.* ,o 1 



E= d 



2^ 2 



-ii 2 
2 



Once again I emphasize that these results are important for the verification of the stability of 
our numerical analysis. This derivation is shown here not only as an exercise, but also because I 
couldn't find the 2-d version in any textbook. 



Chapter 2 

Finite Differences 

2.1 Introduction 

Differently from the non-approximate analytical solutions of PDEs in the continuum (for instance, 
those obtained through variable separation and subsequent integration), numerical solutions ob- 
tained in a computer have limited precision]]. It is due to the way in which computers store data 
and also because of their limited memory. After all, how could we write in decimal notation (or 
in any other base) an irrational number like y/2 making use of a finite number of digits? In this 
work we won't stick with rigorous derivations of the theorems nor of most of the results presented. 
The references listed in the end should be considered for this end. 

The central idea of numerical methods is quite simple: to give finite precision ("the discrete") 
to those objects endowed with infinite precision ("the continuum"). By discretize we understand to 
transform continuum variables like x, y, .., z into a set of discrete values \xi), {j/i}, ..., {zi}, where 
i runs over a finite number of values, thus sampling the wholeness of the original variables. As 
a consequence of this discretization, integrals become sums and derivatives turns out to mere 
differences of finite quantities (hence the name "finite differences"). I illustrate below these ideas: 

I f(x)dx = lim Y^ f{n5x)5x -» ^ f(nAx)Ax (2.1) 

n n 

df(x) = Um f(x + 5x)-f(x) ^ f{x + Ax) - f(x) 
dx 8x^0 5x Ax 

where Sx is a variable with infinite precision (thus its value could be as small as we want) and 
Ax C 1 is a variable with finite precision, which under the computational point of view is the 
limiting case analogous to Sx. We could naively expect that, the smaller the value of Ax, the closer 
we are to the continuum theory. This would be indeed true if computers didn't have finite precision! 
The closer your significant digits get to the limiting precision of the computer, the worse is your 
approximation, because it will introduce the well-known "round-off errors", which are basically 
truncation errors. The reference g] has a somewhat lengthy discussion about computational 
issues like this. 



1 In this point it is worth mentioning that there are basically two ways of solving a mathematical problem 
with the aid of a computer: symbolically and numerically. Symbolic methods deal fundamentally with algebraic 
manipulations and do not involve explicit numerical calculations, giving us an analytical form (whenever possible) 
to the desired problem. It is, however, widely understood that non-linear theories hardly have a closed-form 
solution, and even if they do, it is often a lot complicated and requires an understanding of very sophisticated tools. 
Whenever this is the case, one often resorts to the numerical approach, which doesn't furnish us with an analytical 
closed-form solution, but could give very precise numerical estimates for the solution of the problem. It has been 
used since the very beggining of the computer era, and today it is sometimes the only tool people have to attack 
some problems, pervading its use in almost every discipline of science and technology. 
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2.2 Difference Equations 

Difference equations are to a computer in the same way as differential equations are to a good 
mathematician. That is, if you have a problem in the form of a differential equation, the most 
straightforward way of solving it is to transform your derivatives into differences, so that you finish 
with an algebraic difference equation. This turns out to be necessary for what we said about the 
limitations of a computern. 

As a trivial example, take the ordinary differential equation: 

£=,(,) (2.3) 

Using a first-order Taylor expansion (see Appendix [A|) for f(x), 

f(x + h)nf(x) + f'(x)h^ 

^ f(x + h)-f(x) 

J W ~ h 

we obtain the Euler form for the Eq. ( |2.3| ) : 

f(x + h)-f(x) 
- h 9{x) 

Notice that this equation involves only differences as we said above, and to solve it in a computer 
we shall need the following iterative relation obtained directly from the above equation: 

f(x + h) = hg{x) + f(x) 

or, in the traditional numerical notation: 

In+i = hg n + /„ (2.4) 

Technically, once provided both the initial condition (for instance, /o = 0) and the functional 



form of g„ = g{x n ), we could solve Eq. (2.4) by iterating it in a program loop. 

In spite of its simple form, the Euler method is far from being useful for realistic equations; it 
could give rise to a completely erroneous approximation. Higher order expansions are frequently 
used in order to obtain equations with reduced error (see again Appendix [A] for some of these 
expansions). However, these higher order approximations are also subject to serious problems, 
like the lack of stability or convergence, so the problem is ubiquitous and has been one of the most 
attacked problems in the so-called "numerical analysis", a relatively modern branch of mathematics. 
I will say a little bit more later about these issues on convergence and stability. 

It has already been said that the relevance of the classification of PDEs lies in their "nature"; 
those of initial value have a completely different way of solving numerically from those of boundary 
values. The latter kind doesn't evolve in time. It's the classic case of the Poisson equation which 
could be describing a thermostatic or an electrostatic system: 

d 2 u(x,v) d 2 u{x,y) _ 

dx 2 + dy 2 -H*,V) ( 2 -5) 

Using the expansions from Appendix |A], we have: 

Uj+ij — 2u i: j + Uj-ij Uj,, J+ i - 2ujj + Uij-i 
h 2 h 2 



= /«. (2-6) 



2 There are, however, more sophisticated methods like Finite Elements, but the fact that one needs to get rid 
of differentials transcends these methods when we are talking about numerical solutions. 
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where we took h as the lattice spacing (also grid or net resolution) both for the coordinate x and 
y (x — > x n = nh and y — ► y n = n/i). The indices of this equation then correspond to sites in this 
lattice (a.k.a. lattice points), and sweep from to the number of sites N{ or Nj. The problem 
becomes then to solve the equations given by ( |2.6| ) simultaneously for the various mj. There are 
very interesting methods to solve this sort of problem which could be found in the references || g] . 
Our work, however, is directed towards initial value problems which, as the very name suggests, 
deal with temporal evolutions starting from certain "initial values" at the "zero" instant. It is the 
typical case of the wave equation already presented, or of the diffusion equation 

d 2 u du 

~dtf = W (2 - 7) 

Our focus goes even finer, since we will deal only with "explicit discretizations", which could 
be understood as those which could be solved iteratively, that is, we could solve the difference 
equation for u{x,y,t + At) explicitly in terms of the other variables u(x',y',t) at the instant t 
(for instance, Euler's equation above is an explicit method). Implicit methods need a different 
approach, which often involves the solution of linear systems by using matrices (again [g, || do 
very well in these matters). 

2.3 The von Neumann Stability Analysis 

How should we know, after transforming a differential equation into finite differences, if the calcu- 
lated solution is a stable one? By numerically stable solutions we understand those in which the 
error z 1 ^ between the correct theoretical solution u(x m ,t n ) and the numerical solution U^ does 
not diverge (i.e., is limited) asn^oo (t — > oo), in other words: 

Z rn = llyXm-, tn) U m ^ ^5 V */ 

for any n, where the lower indices are spatial and the upper ones are temporal, and e is a finite 
real value. For instance, an unstable discretization describing a vibrating string could be easily 
detected watching the energy of the system for a while: a divergent energy would certainly arise. 
Fortunately, there is a useful tool to identify unstable finite difference equations prior to simulating 
it, known as von Neumann stability analysis |g, g], which could be applied to a difference equation 
to preview its numerical behavior. 

The von Neumann method consists essentially in expanding the numerical error z^ in a discrete 
harmonic Fourier series: _ 

C = I>(i„)e ifc ^ (2.9) 

r 

and analyzing if a r (t n ) increases (or decreases) as t — * oo (technically, if a r (t) decreases when 
t — ► oo we have a numerical dissipation, which is usually harmless). It is then easy to see that if 
a r (t n ) isn't divergent for any n and m we will have a stable solution. This analysis is somewhat 
simple, since it is sufficient to study the behavior of a single general term of the series, for if we 
prove that this general term of the series could have a certain pathological behavior (like diverging 
for n — > oo), then the whole solution is compromised; otherwise, our solution is stable. 

Mitchell and Griffiths || show that z^ given by ( |2.S| ) satisfy the very same difference equation 
for u^. Hence, if we take a certain z^ such that \z%A = 1 and put it into the difference equation, 
we could achieve the desired stability condition. One possible z^ satisfying the criteria above is: 

C = e anAt e ll3mAx (2.10) 

Indeed, notice that for n — we have |z^| = 1, and with a and arbitrary values we satisfy 
the above discussion. With this expression, we could now write the stability condition for the von 
Neumann analysis: 

\C\ < 1, (2.11) 
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12 



Figure 2.1: Representing an integration algorithm that needs the values u(x — Dx,t) and u(x + 
Dx,t) to obtain u(x,t + Dt). The ratio Dx/Dt is then the "maximum speed" with which an 
information in the algorithm could propagate. 



where £ = e™ is the amplification factor. In summary, putting the error given by 

n fn if3mAx 



(2.12) 



into the difference equation, together with (2.11), we get the necessary condition for stability. 



2.4 The Courant Condition 



Another important condition that we should pay some attention in initial value problems is related 
to the speed with which information could propagate in the difference equation. We could visualize 
the problem in the scheme of Figure 2.1. 

It could be shown || that, applying von Neumann's condition for hyperbolic problems we 
arrive at the Courant condition: if the "physical" wave velocity \v\ in a differential equation is 
greater than the "algorithm speed" Ax/ At, then the scheme is unstable. We have therefore the 
following expression for the Courant condition: 



H< 



Ax 
At 



(2.13) 



2.5 The Staggered Leapfrog Algorithm 

In this section we present the basic idea behind the staggered leapfrog method, which is a very 
efficient second-order integration scheme without numerical dissipation. 

The method is for differential equations which could be cast on the form of a flux-conservative 
one: 



(2.14) 



du _ dF 

dt dx 

For instance, the wave equation could be written in this form when we make the following substi- 
tutions (to ease the notation, from now on we will take v = t = A = a = 1): 



d du 
dt'dt ~~ 


d du 
dx dx 


ds dr 


ds dr 


dt dx' 


dx dt 



where s = dtu and r = d x u. So, putting in the required form, we have: 



d_ 
dt 



d_ 

dx 



— r 
—s 
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The name "leapfrog" refers to the fact that both the time and spatial derivatives are to be 
evaluated using a centered difference scheme, so the integration is done from time-step n — 1 to 
n + 1, "leaping" over the spatial derivatives which are evaluated at step n. The centered difference 
is a second-order one, giving for the flux-conservative equation: 

U i U i _ r i+l r i-l 

2At 2Ax 

In the case of second-order differential equations, we need also the first derivative in order 
to integrate for the next step. In the present scheme, these derivatives are evaluated at some 
"artificial" points between the "real" ones, and hence are said to be "staggered" with regard to the 
variable u n which is evaluated at the real points. Take for example the case below in which an 
arbitrary function f[u] does not depend on any time derivative of u: 

d -^- = fHx,t% (2.15) 

where s = dtu and, for the wave equation, f(u) = d 2 u. According to the above discussion, the 
derivative s is to be evaluated at points n + 1/2, n + 3/2, ..., so we haveQ 



-1/2 _ 9u 
dt 



1+1/2 u n+l 



At 



u «+! = u » + Ats™ +1/2 (2.16) 



But Eq. ( 2.15 ) furnishes us with an expression for s 

n+l/2 n-l/2 



n+l/2_ 



At 



= f{u n ) 



s r +l/2 = *r 1/a + At /K i ) 

We could now write the staggered leapfrog method for equations like ( [2.15] ): 

s „+l/2 = s n-l/2 + Atf{un) (^V 

It could be easily shown Q that this procedure is formally equivalent to taking direct second- 
order finite differences for the second-order derivatives in the wave equation, giving the following 
scheme: 

u i+l ^u t i-%_i _ u % zu i + u i 

Ax 2 At 2 \ ■ ) 

The situation with a dissipative term is a little bit more complicated. We said above that the 
first derivatives in the staggered leapfrog method are to be evaluated at the staggered points, but 
an equation with a dissipative term has a first derivative evaluated at the real ones: 

^ = -VS + f[u(x,t)] ^ (2.19) 

n+l/2 _ n-l/2 

Sj At St = ~VS? + /[<] 



Notice that the centered differences are now evaluated with a step of At/2. 
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Indeed, the above equation has an explicit first derivative evaluated at a real point (sf). How 
can we fix this problem? The trick is to take an average over the adjacent points, so we recover 
the staggered points in the dissipative term: 



n+1/2 n— 1/2 n+1/2 . re-1/2 

Sa — s, s„- + s.- 



Af 



7K 



Continuing with the same procedure done for the previous case (we need an expression for 

1+1/2 



to use with Eq. (2.16)), we have: 



, + 1/2 + ^ s „ +1/2 = s „_ 1/2 _ ^n-1/2 + Ai/[<] ^ 



„ + i/ 2 (i-^)»r 1/a +**/[* 

S, = 



2 



Now we have the staggered leapfrog scheme for equations like ([2.19|): 



•s 



u? +l = u? + Ats7 +1/2 



1 + ^ 



For the case of a wave equation in d dimensions, the f[u] term will have a d-dimensional 
Laplacian V 2 u(x,t), with x — (x\,xi, ...,Xd), and this Laplacian is to be evaluated numerically 
with a second-order finite difference scheme at the step n: 

„2 /-> .^ a i+l,j.k,... za i,j.k.... + a i-l,j,k,... , u i,j + l.k,... ^ a i.j,k,... ' "i.j-l.fc,... . 
VU{X,t)Ki ' 



Axi A^2 

1/2 

When solving the staggered leapfrog method for the first step n = 1, we need s i / , but the 
initial conditions are usually defined at the initial time n = 1. One way to solve this problem is 
to integrate the first step using the Euler scheme (which needs only the initial points at n = 1) 

■I /o 

with a time-step of At/2, thus obtaining the required value of s i . 

Now some words about the efficiency of this method. When we apply the von Neumann 
stability analysis in leapfrog equations we arrive at the Courant condition || . When this condition 
is satisfied, the staggered leapfrog method is not only stable, but also conservative, in the sense 
that it does not introduce any numerical dissipation. 

For the applications which will be shown in Chapter |l we shall use the following conditions: 



du(t) 
u ^' I boundaries — U ' ~~Qt~ 



= 
boundaries 



(r — rj)) 2 
u{x,y)\ t=0 = Cexp , 

where f = xi + yj, r*o = hi + |j, I is the lattice length, C is a normalization constant, and 7 is a 
sufficiently small constant such that u(r) — * as r -^boundaries, that is, the initial condition is a 
gaussian sufficiently localized to make u smooth at the boundaries. 



Chapter 3 

Examples 



3.1 The Free String (ID) 

These simulations were executed in a PC of 350MHz, for lattices of at most N = 1000. The 
integration time lay in the order of seconds. 



Figure 3.1 shows some results for the conditions of the previous section for various parameters. 
Notice the improvement of energy conservation for finer resolutions and the very good exponential 
fitting for i] — 1. This exponential result is expected, since the vibrating string could be understood 
in terms of the Fourier space, where each mode behaves as a decoupled harmonic oscillator with 
a damping given by r\. 



3.2 The Membrane (2D) 



These simulations were also executed in a PC of 350MHz, and for lattices of N = 500 the 
integratio n tim e reached half an hour. 



Figure |3,2| shows some results for r\ — and r\ = 1. For a non-conservative system (r\ — 1), we 



see also the exponential fitting for N = 200. 
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Figure 3.1: E x t for different lattice spacings and r\ = (above) and E x t for N = 1000 and 
r] = 1 (below). The linear dimension is fixed at L = 1. 
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Figure 3.2: E x t ior ij — and various TV (above) and E x £ for r\ — 1 and JV = 200 (below). The 
dimensions of the membrane are L x L = 1 for all runs. 



Appendix A 

Taylor's Theorem 



A.l Definitions 

When we want to transform a differential equation into a difference equation, the Taylor expansion 
is often used: 

fix) = f(x t) ) + f(xo)(x-xo) + ^f"(x )(x-x Q ) 2 + ...= (A.l) 

. tqpl {x _„,., (A „ 

n=0 

where so is the point around which we want to expand /(£)[]. An alternative form for this 
expansion is achieved doing a simple variable change x —* x + h and xq — > x: 

00 £ In) ( \ 

f( x + h) = J2 L ^f l h n (A-3) 

* — ' n! 

In the numerical case we will be interested in truncating the series, so that we finish with a 
finite number of terms. We could then write this expansion in the following form: 

m f(n) I \ 

f( X+ h) = y2^^h n + o(h m+i ), (a.4) 

n=0 

where 0(h m+1 ) corresponds to the truncated terms which powers of h are equal or higher than 
TO+1 (this term is frequently called the error of order m + 1). Notice that under the numerical 
point of view it is important to know the order of O in the discretization, since for h -C 1, the 
greater the order of O the more negligible will be the error. 

A. 2 Useful Expansions 

Some expansions that will be used throughout this text are shown below. All of them could be 



obtained from (A.4) by directly solving for the desired term or using more than one expansion to 



find higher order expansions for the derivative, and then solving the system. For instance: 

fix + h) = f(x) + f(x)h + y"{x)h 2 + ... 

fix - h) = f{x) - f'(x)h + \f"(x)h 2 - ... 

fix) = f{X + h) - f{X) -0{h) (A.5) 



1 For xo = this expansion is also known as Maclaurin expansion. 
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fix) = /(* + *) -A*-*) - 20(ft 2 ) (A.6) 

/»(*) = /( - + fe) - 2/ ^ ) + /(X ' fe) - 20(*») (A.7) 



df{x,y) _ f{x + h,y)- f{x 2 y)__ 



dx h 

d 2 f{x 1 y) _ f(x + h,y)- 2f(x, y) + /(.x - ft, y) 
9a; 2 ft 2 



2C(ft 2 ) (A.9) 



Notice that when we divide an error of order 0{h n ) by h r , automatically this error turns to order 
n - r, i.e., 0{h n )/h r = 0(h n - r ). 
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